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Technical  Management  Office  (TMO),  Science  &  Technology  Group  (S&TG), 
Computational  Analysis  Team  (CAT). 

The  author  wishes  to  thank  Mr.  Steven  E.  Kunz  of  S&TG's  Experimental 
Research  Team  (ERT),  for  his  helpful  suggestions  concerning  the  thermodynamics 
of  the  retraction  process,  and  providing  the  Tustin  experimental  data  used  in 
comparing  model  predictions  with  actual  system  behavior. 


PREDICTIVE  MODEL  OF  A  PARACHUTE  RETRACTION  SOFT 

LANDING  SYSTEM 


Introduction 

The  U.S.  Army  Soldier  Systems  Command's  Natick  Research,  Development 
and  Engineering  Center  (NRDEC)  is  currently  examining  a  concept  for  reducing 
the  impact  shock  sustained  by  airdropped  payloads  upon  ground  impact.  A  device, 
called  a  parachute  retractor,  is  placed  between  the  payload  and  parachute 
confluence  point,  and  when  activated,  accelerates  the  parachute  and  payload  toward 
each  other;  slowing  the  payload  prior  to  ground  impact.  The  goal  is  to  eliminate 
the  cushioning  material  currentiy  placed  under  airdrop  loads,  and  thus  eliminate  the 
time-consuming,  labor-intensive  process  of  rigging  payloads  for  airdrop,  while 
providing  a  roll-on/roll-off  (RO/RO)  capability  for  vehicles. 

The  concept  of  parachute  retraction  to  soft  land  cargo  is  not  new.  Its  origins 
date  back  to  a  1956  report  on  airdrop  cushioning  prepared  by  the  University  of 
Texas  for  NRDEC,  then  known  as  the  Quartermaster  Research  and  Development 
Command^  In  that  report,  the  characteristics  of  an  elastic  spring-type  device  for 
use  with  parachutes,  dubbed  a  landing  snubber,  were  described  and  a  theoretical 
analysis  of  its  operation  presented.  A  further  expansion  of  the  concept,  known  as 
Parachute  Reel-In  Reel-Out,  was  examined  by  NRDEC  in  the  1960's.  In  addition 
to  soft  landing  cargo,  this  concept  envisioned  use  of  retraction/urureehng  of  the 
hues  connecting  die  parachutes  and  payload  in  order  to  maintain  high  relative 
velocities  of  the  airstream  with  respect  to  the  parachutes,  and  to  relieve  unwanted 
high  forces  on  the  cargo  at  appropriate  times  during  airdrop  descent.  Use  of  a 
powered  winching  device  located  between  the  recovery  parachutes  and  the  cargo 
to  eitiier  decrease  or  increase  the  distance  between  the  parachutes  and  cargo  during 
the  trajectory  was  envisioned.  None  of  these  concepts  ever  reached  the  prototype 
stage  or  were  tested. 

In  the  spring  of  1994,  instrumented  airdrop  tests  of  a  novel  retraction  device 
invented  by  John  Lanza  of  Natick  were  conducted  for  the  first  time  at  Tustin 
Marine  Corps  Air  Station  (MCAS),  Santa  Ana,  CA.  Experimental  data  obtained 
in  these  tests  showed  that  the  Natick  parachute  retractor  could  slow  or  even  stop 
payloads  just  before  landing.  The  retractor  consisted  of  a  pneumatically  driven 
piston  connected  by  cables  to  upper  and  lower  pulley  blocks;  all  enclosed  within 
a  cylinder.  The  pulleys  are  used  to  increase  the  length  of  retraction  compared  with 
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the  stroke  of  the  piston,  as  is  shown  schematically  in  Figure  1. 

The  benefit  of  developing  a  model,  which  can  predict  the  performance  of  this 
Parachute  Retraction  Soft  Landing  System  (PRSLS)  for  all  airdrop  configurations, 
including  large,  yet  to  be  investigated  (experimental)  systems,  was  recognized  early 
in  the  program.  Therefore,  as  part  of  this  exploratory  development  effort,  a 
computational  model  was  developed  to  determine  the  optimal  activation  height,  to 
conduct  system  trade-off  studies,  and  to  investigate  scaling  effects.  This  report 
reviews  the  development  of  that  computer  model. 

The  report  begins  with  a  review  of  the  model's  underlying  theory  and 
assumptions.  Applying  the  principles  of  thermodynamics,  kinetics,  kinematics, 
fluid  mechanics  and  numerical  methods,  key  mathematical  relationships  used  in 
modeling  the  system  as  a  three-mass,  three-degree-of-ft-eedom  system  are  first 
derived.  Derivations  are  given  exhaustive  treatment  in  order  to  enhance  the 
reader's  understanding  of  the  theory,  principles,  and  assumptions  that  provide  the 
foundation  for  model  development.  Rigorous  treatment  is  given  to  the  development 
of  expressions  for  the  forcing  funtions,  i.e.,  force  acting  on  the  piston,  the  work 
done  by  the  expanding  gas  on  the  piston,  the  differential  equations  of  motion  of  the 
discrete  masses,  the  governing  nonlinear  differential  equation  of  motion  of  the 
system,  and  the  recurrence  formulas  used  to  integrate  the  motion  equation. 
Solution  of  the  equation  of  motion  and  the  trial  and  error  iterative  scheme  used  to 
reconcile  the  system’s  thermodynamics  and  kinetics  and  arrive  at  the  time  values 
are  discussed  next.  Finally,  FORTRAN  computer  model  predictions  are  compared 
to  die  data  obtained  during  experiments  conducted  at  Tustin  MCAS  on  20  April 
1994,  and  conclusions  are  drawn  concerning  model  accuracy. 
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To  Parachute 
Confluence  Point 


Parachute 


Thermodynamics 


The  principles  of  thermodynamics^’^  are  applied  in  this  section  to  the 
development  of  expressions  for  the  force  exerted  on  the  piston,  and  the  work 
produced  by  the  expanding  gas,  upon  valve  opening.  The  piston/cylinder 
mechanism  is  assumed  to  be  frictionless.  Force  and  work  expressions  are  derived 
as  a  function  of  piston  displacement,  y,  for  the  processes  of  both  adiabatic  and 
isothermal  ideal  gas  expansion  within  the  cylinder.  Constants  include  the  tank 
volume,  Vj  ,  initial  tank  pressure,  Pj  ,  specific  heat  ratio,  k,  (adiabatic  only)  and 
piston  radius,  r.  Derivations  for  each  process  begin  with  an  examination  of  die 
First  Law  of  Thermodynamics,  which  in  differential  form,  is  written  as  follows: 

dQ=dW^dU 


This  states  that  the  heat  entering  die  system,  dQ,  is  equivalent  to  the  work  done 
against  the  atmosphere,  dW,  plus  the  internal  energy  accumulation,  dU. 

Adiabatic  Expansion  Process 

An  ideal  gas  expands  from  initial  state  (1)  to  final  state  (2)  in  a  cylinder 
which  is  completely  insulated  so  that  there  is  no  heat  exchange  with  the 
surroundings;  that  is:  Q=0. 

The  First  Law  of  Thermodynamics  in  differential  form  is  written  as: 

dQ=dW+dU=0 


or 


dW^=-dU 

which  indicates  that  the  work  effect  is  produced  at  the  expense  of  the  system's 
internal  energy. 

For  an  ideal  gas,  the  enthalpy,  H,  which  is  equal  to  U+PV,  and  internal  energy,  U, 
depend  on  the  temperature  only  and  are  independent  of  pressure. 
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Therefore,  the  following  relationships  hold  true  at  all  pressures: 


where  Cp  and  are  respectively  the  specific  heats  at  constant  pressure  and  constant 
volume,  and  R,  the  universal  gas  constant,  is  equal  to  Cp-c^.  It  should  be  apparent 
from  the  above  relationships  that  for  isothermal  processes  (i.e.,  constant 
temperature,  AT=0),  both  the  change  in  enthalpy,  AH,  and  internal  energy  change, 
AU,  are  equal  to  zero.  Returning  to  the  adiabatic  case,  it  follows  that  for  gas 
expansion  within  the  cylinder: 


or 


dW,^PdV=-nc,<iT 

where  n  is  the  number  of  moles  of  gas  present.  Upon  differentiating  the  Ideal  Gas 
Law  equation: 

PV=nRr 


one  obtains  the  following: 

PdV+VdP=nRdT 


Upon  substituting  (-PdV/nCy)  for  dT  in  the  above  equation  one  obtains: 

PdV+VdP=-nRi 


^PdV^ 


Rearranging  terms  results  in  the  following  expression: 
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But  R=Cp-Cy.  Therefore: 

c  V  P 
\  vj 

Integrating  between  states  (1)  and  (2): 


where  k  =  Cp  /  c^. 


•••  PiVi^p^Vi 


It  should  be  noted  that  for  a  Monatomic  Gas: 

Cp  =  5/2  R 
c,  =  3/2  R 

and  for  a  Diatomic  Gas: 


Cp  =  7/2  R 
c,  =  5/2  R 


Since  PjVi  =  nRTi  and  P2V2  =  nRT2,  upon  substitution  of  nRTj  /  Vj  and 
nRT2  /  V2  respectively  for  Pj  and  P2  in  the  above  equation  one  obtains: 


(nRT> 

II 

rnsr^'i 

i  J 

‘  1 

(1) 


Similarly,  if  one  substitutes  nRT,  /  Pj  and  nRT2  /  P2  respectively  for  Vj  and 


N 2  01^^  obtains: 


(nRT^ 


\k 


\  ) 


=P. 


nRT^ 


y 


Px  P2 


1 


pt-1  pt~l 
M  ^2 


/T 


Pi 


\k-l 


T 

.  2 
■  Pi 

The  adiabatic  Volume-Temperature  and  Pressure-Temperature  relationships 
shown  respectively  in  (1)  and  (2)  were  not  incorporated  into  the  model.  These 
equations  were  derived  and  are  presented  nonetheless,  since  they  would  be  needed 
to  develop  expressions  to  predict  the  gas  temperature  from  valve  opening  up  to  the 
full  piston  stroke,  if  desired  at  some  later  time. 


Given  an  initial  tank  pressure,  Pj  ,  and  tank  volume,  ,  the  gas  volume  at  some 
later  time  during  the  expansion  is  given  by  the  expression: 

K(y)=P)+7cr^ 


where  r  is  the  Piston  Radius  and  y  is  the  Piston  Displacement. 

From  the  previously  derived  Pressure-Volume  relationship  for  adiabatic  ideal  gas 
expansion  it  is  evident  that: 


Upon  substitution  of  the  above  expression  for  V(y)  one  obtains  the  following: 

PV^ 

P(y)=-  ‘  ‘ 


Noting  that  the  force,  F,  acting  on  the  piston  is  equal  to  7cr^P(y)  it  follows  that: 
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(3) 


nr^Pyi 

F(y)=- 


But  dV/dy  =  7cr^  or  dV  =  Tcr^dy.  Therefore,  the  work  done  on  the  piston  during  the 
adiabatic  expansion  of  the  gas  within  the  cylinder  is  given  by  the  following: 


rs  ,  rs  PiVi 

W=f  PdV=f  - 5-^— - 

«  Jo  Jo 


where  S  is  the  piston  stroke  and  0  <  y  <  S. 

Upon  collecting  and  rearranging  terms,  one  obtains  the  following  expression: 

W=Kr'^PyH^ - ? - dy 


or 


W^=-nr^P,V^fJ{V,*nr^}-^dy 


Let  u  =  Vj  +  Tcr^  .  Then  du  =  Tcr^dy  or  dy  =  du  /  Tcr^ .  Upon  substitution  into  the 
above  equation  for  one  obtains  the  following: 


Integrating  the  above  expression  results  in: 


W  =PV 

q  *  r  I 


u 


1-jt 


1-^ 


S 

0 


Substituting  +  Tcr^y  for  u  in  the  above  equation  results  in: 
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s 


W  =P.V 

a  ■*  r  i 


(1-t) 


Evaluating  the  above  integral  results  in  the  following  equation: 


W  =P.V 

”q  y  I 


(P,V,] 


1-k 


1-k 


PiV,\Vi+nr^Sf-^-P.Vi 


The  total  adiabatic  work  done  by  the  expanding  gas  up  to  some  displacement,  y, 
during  piston  movement,  is  similarly  obtained  by  evaluating  the  above  integral 
between  zero  and  the  value  of  y.  Thus: 


PiV-(V,^nr^f-^-P,V, 

(l-k) 


(4) 


where  0  <  y  <  S. 

Isothermal  Expansion  Process 

An  ideal  gas  expands  from  initial  state  (1)  to  final  state  (2)  in  a  cylinder. 
The  temperature  remains  constant  throughout  the  process.  The  volume  of  the  gas 
increases  as  the  pressure  of  the  gas  decreases.  The  cyhnder  is  in  contact  with  a 
constant  temperature  bath  which  adds  heat  to  the  system  as  the  system  does 
expansion  work. 

The  First  Law  of  Thermodynamics  requires  that: 

dQj,=dW^dU 


Since  temperature  is  constant,  i.e.,  dT=0,  dU  of  the  ideal  gas  is  equal  to  zero,  and 
all  the  heat  entering  the  system  leaves  as  work.  Therefore: 
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dQ^^dW 


and  the  Ideal  Gas  Law  in  differential  form  is: 

PdV^VdP=Q 


Therefore: 


V  P 


Integrating  between  states  (1)  and  (2); 


inKi;;;=-inPi;:; 

(lnK,-lnK,)=lnF,-lnPj 


(Vo] 

In 

2 

=ln 

1 

1  n 

ii-£i 

K.  P, 

••• 

Given  an  initial  tank  pressure,  ,  and  tank  volume,  Vj ,  the  gas  volume  at  some 
later  time  during  the  expansion  is  given  by  the  expression: 

V(y)=Vi+nr^ 

where  r  is  the  Piston  Radius  and  y  is  the  Piston  Displacement. 
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From  the  Ideal  Gas  Law  for  a  constant  temperature  process,  it  is  evident  that: 


Substitution  of  the  previously  derived  expression  for  V(y)  into  the  isothermal  Ideal 
Gas  relationship  results  in  the  following: 


PV. 

P(y)=-  ‘  * 


Noting  that  the  force,  F,  acting  on  the  piston  is  equal  to  Ttr^PCy)  it  follows  that: 

nr^Py, 


I  I 


(V^+Tir^) 


(5) 


From  die  volume  expression  it  follows  that  dV/dy  =7cr^  or  dV  =  jcr^dy.  Therefore, 
the  work  done  on  the  piston  during  the  isothermal  expansion  of  the  gas  within  the 
cylinder  is  given  by  the  following: 


where  S  is  the  Piston  Stroke  and  0  <  y  <  S. 

Let  u  =  Vi  +  TcrV-  Then  du  =7a^dy  or  dy  =  du  /  Tcr^.  Upon  substitution  into  the 
above  equation  for  one  obtains  the  following: 

W^=P,V,KrH^^-^=P,V^  f 

•'«  «  icr^  ' 

Substituting  the  expression  for  u  back  into  the  above  equation  and  evaluating  the 
integral  results  in  the  following: 
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w^=py^ 


(  V,+nr^S 


The  total  isothermal  work  done  by  the  expanding  gas  up  to  some  displacement,  y, 
during  piston  movement,  is  similarly  obtained  by  evaluating  the  above  integral 
between  zero  and  the  value  of  y.  Thus: 


/ 

W^)=:P,V^  ■ 
\ 


V^+nr^ 


(6) 


The  above  work  expression  can  be  derived  alternatively  by  noting  that  the  Ideal 
Gas  expression  PV  =  nRT  permits  the  elimination  of  either  P  or  dV  from  the  work 
expression: 

W=f^PdV 


From  the  Ideal  Gas  Equation  it  is  apparent  that: 


nRT_RT 
V  ■  V 


or 


where 


dV=—dP+-dT 

p2  p 


n 

is  defined  as  the  volume  per  mole  or  molar  volume. 
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Upon  substitution  of  the  above  expressions  into  the  work  equation,  one 
obtains  the  following: 


or 


But  dT  =  0.  Therefore: 


^dP+MT^ 


Let  PjVj  =  nRT  =  PjVi.  One  will  note  that  the  judicious  substitution  of  Pj ,  V;  and 
the  previously  derived  expressions  for  V(y)  and  P(y)  respectively  for  Pj  ,  Vj  , 
and  P2  in  either  of  the  above  work  expressions  results  in  the  previously  derived 
work  equation  (4)  for  isothermal  gas  expansion  in  a  cylinder. 
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Development  of  the  Differential  Equations  of  Motion 

The  PRSLS  was  modeled  as  a  three  mass,  three  degiee-of-freedom  (DOF) 
rigid  body  system.  Assumptions  made  in  order  to  simplify  modeling  of  the  system 
included  the  following: 

(1)  Piston  and  pulley  movement  is  frictionless. 

(2)  Pulley  cables  are  inextensible,  i.e.,  cable  length  remains  constant. 

(3)  System  orientation  remains  vertical  throughout  the  retraction. 

(4)  The  apparent  mass  of  the  parachute,  which  consists  of  air  inside  the  canopy, 
as  well  as  that  in  the  immediate  vicinity,  is  constant  and  equal  to  the  mass  of  air 
contained  within  a  hemisphere  of  nominal  parachute  diameter. 

(5)  The  parachute  drag  coefficient,  based  on  the  nominal  parachute  diameter, 
remains  constant  during  retraction. 

The  PRSLS  model,  as  previously  noted,  is  comprised  of  three  discrete 
masses.  They  include:  (1)  the  piston  and  any  upper  retraction  device  cross 
members  and  upper  pulley  blocks  that  move  along  with  the  piston;  defined  as 
MpisT  ,  (2)  the  payload  and  the  retraction  device  cylinder,  tank,  lower  pulley  blocks 
and  structure  that  move  with  the  payload;  defined  as  Mp  ,  and  (3)  the  parachute 
with  mass  concentrated  at  the  parachute  confluence  point;  defined  as  M^.  The 
three  degrees-of-ffeedom  (DOF)  considered  are  the  vertical,  i.e.,  Y-direction, 
displacements  of  the  piston,  payload  and  parachute  confluence  point;  defined 
respectively  as  y,  yp  and  y^.. 

The  retraction  device  modeled  contains  two  banks  of  upper  and  lower 
pulleys,  and  two  sets  of  pulley  cables.  Each  of  the  pulley  cables,  which  converge 
at  the  parachute  confluence  point,  by  symmetry  carry  one-half  the  tension,  T^,  at 
the  parachute  confluence  point  at  any  instant  in  time  during  the  system's  descent, 
just  prior  to  and  during  retraction.  Formulation  of  the  differential  equations  of 
motion  was  carried  out  by  applying  the  concept  of  dynamic  equilibrium  or 
Newton's  second  law  of  motion  to  the  system"^.  For  tiiis  purpose,  the  three  masses 
were  first  isolated  showing  all  of  the  forces  acting  on  each  mass,  including  internal 
forces,  externally  apphed  forces  and  inertia  forces.  A  detailed  examination  of  the 
forces  acting  on  fire  piston  and  payload  masses,  i.e.,  Mpi^  and  Mp  ,  revealed  two 
important  factors  concerning  internal  forces  developed  within  the  retractor's  pulley 
subsystem.  First,  for  a  pulley  arrangement  of  a  given  mechanical  advantage,  N,  the 
force  on  the  piston  tending  to  resist  piston  movement  is  equal  to  2  x  N  x  (JJ2) 
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=  N  X  Tc-  Second,  for  a  mechanical  advantage,  N,  the  force  tending  to  accelerate 
the  payload  upward  is  equal  to  2  x  (N+1)  x  (JJ2)  =  (N+1)  x  Tc- 


Applying  Newton's  second  law  of  motion,  the  equation  of  equilibrium  for  the 
piston  in  the  vertical  direction,  +t  E  =  ma,  is  written  as: 

(7a) 

Similarly,  the  equihbrium  equation  for  the  payload  in  the  vertical  direction, 
+TEf,  =  ma,  is  written  as: 

{N^\)Tc-F-Wp-Mp%  Ob) 


The  force,  F,  in  the  above  equations  of  motion,  varies  with  time,  i.e.,  F  =  F(t),  and 
is  the  force  applied  to  the  piston  by  the  expanding  gas  upon  valve  opening. 

Finally,  applying  the  principles  of  dynamic  equilibrium  at  the  parachute 
confluence  point  in  the  vertical  direction,  +i  E  F^  =  ma,  results  in  the  following 
equilibrium  equation: 

Oc) 


The  parachute  drag  force  and  apparent  mass  temnis,  and  Me  ,  in  the 

above  equation  are  given  by  the  expressions^: 


TirpCjjP 


Parachute^ 


<?cf 


(8) 


and 


Mc=^T^prl  (9) 


where  rp  is  the  nominal  parachute  radius,  which  is  equal  to  one-half  the  nominal 
parachute  diameter,  Cp  is  the  drag  coefficient  based  on  the  total  canopy  surface 
area,  and  p  is  the  density  of  air.  In  order  to  simplify  analysis,  the  parachute  model 
therefore  assumes  that  both  the  apparent  mass  and  chag  coefficient  remain  constant 
during  retraction. 
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If  one  defines  the  constant,  Cj  ,  where: 


C  = 


nrpC^p 


(10) 


Then  it  follows  that: 


(11) 


Substituting  the  above  expression  for  Fparachute  (7c)  and  noting  the  weight 
of  the  air  in  the  parachute,  ,  is  negligible  and  can  be  ignored,  one  obtains  the 
following: 

(12) 

Substituting  the  above  expression  for  Tc  into  (7a)  and  (7b)  results  in  the 
following  expressions: 

(13o) 

(136) 

Dividing  both  sides  of  (13a)  and  (13b)  respectively  by  Mpjgj  and  Mp  result 
in  the  following  expressions  for  the  piston  and  payload  accelerations: 


,_\F-NC^i^cf-NM^c-^PisT 


M, 


(14a) 


PIST 


yp 


M. 


(Ub) 


The  next  step  is  to  define  the  retraction,  retraction  rate,  and  retraction 
acceleration  in  terms  of  the  displacements,  velocities,  and  accelerations  at  the 
piston,  payload  and  parachute  confluence  point.  Based  on  the  sign  convention  used 
in  deriving  the  dififerential  equations  of  motion,  the  following  relationships  are 
defined: 
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L=Ny=yc+yp 

(15a) 

L=Ny=yc-^yp 

(15i») 

L=Ny=yc^yp 

(15c) 

Therefore,  it  is  evident  that: 

yp=^-yc 


Substituting  the  above  expression  for  the  payload  acceleration  into  (13b)  and 
dividing  both  sides  of  (13a)  and  (13b)  respectively  by  Mpjsx  and  NMp  results  in  the 
following  expressions  upon  collecting  and  rearranging  terms: 


iVCi 


Mpigr  Mpjgj. 


i^cf 


NM. 


M, 


-yc 


PIST 


pm 


M, 


=3^ 


PJST 


(16a) 


NMp 


(iV+l)M.+Mp 


NM, 


F 

NMp 


^-y 

NMp 


(m 


Equating  (16a)  and  (16b)  results  in  the  following  expression: 

(N+l)Mc+Mp 


NC 


NM. 


Mpj^  Mpjgj. 


1  /.  ^2_  ..  _  _(^'^l)^l,.  ,2 


pm 


NM. 


NM. 


P^c- 


F  g 

NMp  N 


where  Mpisxg  and  Mpg  have  been  substituted  respectively  for  Wpigj  and  Wp, 
cancelling  out  the  mass  term  in  the  last  entry  on  both  sides  of  the  expression. 


The  following  nonlinear  differential  equation  of  motion  for  the  parachute 
retraction  soft  landing  system  is  obtained  upon  collecting  and  rearranging  terms  in 
the  previous  expression: 
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[(iV+DCi  JVC,],.  JVJlfcl..  (  g\  J  1  1  ] 

[  NM^  [  JVJlf,  N)  (J\f,,s,.  JVJlf,j 

Multiplying  both  sides  of  the  equation  shown  above  by  NMpMpjgj  results  in 
tile  following  motion  equation: 

Define  the  following  constants: 

C^={m\)C,M^jsT^N^C,M^ 

(18) 

=NMpMpjg-^  -MpMpjgjg  ={N- 1  )MpMpjgjg 

Cj  =NMp  '*'Mpjgj 

Then  the  nonlinear  differential  equation  of  motion  shown  in  (17)  can  be 
expressed  as  follows: 
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Linear- Acceleration  Method 


Integration  of  the  nonlinear  differential  equation  of  motion  developed  in  the 
previous  section  was  carried  out  using  the  Linear-Acceleration  method^’^.  As 
implied  by  its  name,  this  method  assumes  that  the  acceleration  varies  linearly 
between  time  stations  as  shown  in  Figure  2. 


acceleration,  x 


Figure  2:  Linear-Acceleration  Method;  Actual  vs.  Assumed  Acceleration 


The  acceleration  between  time  stations  s  and  s+1  would  therefore  be 
approximated  by  the  expression: 


At 


where  the  superscript  enclosed  by  parentheses,  by  definition,  represents  the 
particular  time  station. 
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The  velocity  at  any  time  within  this  interval  may  be  obtained  from  the 
following: 


/,> 


2At 


which  at  station  s+1  becomes: 


^(s+l)=^(s)  +^^-(s+l)+jf  (s)) 


The  displacement  at  s+1  is  given  by: 


J,(S) 


6 


where  from  the  general  velocity  expression: 


The  recurrence  formula  shown  above  can  be  expressed  as  follows: 

(21) 


where 
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(22) 


2 


From  the  previous  discussion  it  can  similarly  be  reasoned  that: 

6 

This  equation  may  be  expressed  in  terms  of  the  acceleration  at  time  station  s  as 
follows: 


X(s) 


(Atf  At 


Substituting  the  above  expression  for  the  acceleration  at  time  station  s  into 
the  previous  equation  for  the  velocity  at  time  station  s  results  in  the  following: 

At  2 


or 


At  At  2 


and 


At  At  2 


Therefore 


The  above  recurrence  formula  can  be  expressed  as: 
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where 


(24) 

3 


^(5-1)  =_;p(s-i)+2je(«-i)^  +jg(^-i)i^  (25) 

3  6 

In  applying  this  method,  equations  (21)  and  (24)  are  used  in  succession. 
Since  the  vdue  of  the  acceleration  at  time  station  s  is  not  known  in  advance,  the 
solution  must  be  iterative  within  each  step.  This  iterative  process  is  not  self¬ 
starting,  however,  since  it  requires  a  supplementary  formula  to  determine  the  first 
estimate  of  the  velocity  at  time  station  s  in  each  time  step.  First  estimates  of  this 
velocity  were  obtained  using  the  formula: 

(26) 

in  which  die  acceleration  is  assumed  constant  during  the  time  interval  and  equal  to 
the  initial  value. 
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Integration  and  Solution  of  the  Nonlinear  Differential 
Equation  of  Motion 

Using  the  thermodynamic  and  kinetic  relationships  and  numerical  recurrence 
formulas  derived  in  the  previous  sections,  a  FORTRAN  computer  model  was 
developed  to  perform  a  time  history  analysis  of  the  Parachute  Retraction  Soft 
Landing  System  (PRSLS).  This  model  uses  a  two-step  iterative  procedure  to 
integrate  the  nonlinear  differential  equation  of  motion  in  time;  providing  a  time 
history  of  the  motion  of  the  piston,  payload  and  parachute  confluence  point,  i.e., 
displacements,  velocities  and  accelerations,  force  exerted  on  the  piston  by  the 
expanding  gas,  gas  expansion  work,  line  tension,  parachute  force,  and  retraction, 
retraction  rate  and  retraction  acceleration  from  valve  opening  up  through  full  piston 
stroke.  The  FORTRAN  coding  for  this  computer  model  is  listed  in  Appendix  A 
of  this  report. 

Model  inputs  include  the  initial  tank  pressure,  tank  volume,  piston  radius, 
specific  heat  ratio,  weight  of  piston  including  cross  member  and  top  pulleys, 
payload  weight  including  retractor,  piston  stroke,  step  increment  of  piston 
displacement  at  which  time  history  values  are  desired,  gas  expansion  process,  i.e., 
adiabatic  or  isothermal,  parachute  diameter,  mechanical  advantage  of  pulley  system, 
and  parachute  drag  coefficient. 

Using  the  above  model  inputs,  the  program  begins  the  analysis  by  first 
computing  the  three  mass  terms,  the  terminal  descent  velocity  of  the  system,  and 
the  constant  terms  used  in  the  differential  equation  of  motion.  The  mass  at  the 
parachute  confluence  point  is  computed  using  equation  (9).  Values  of  the  constant 
terms.  Cl  through  C5,  are  computed  using  equation  (10)  and  equations  (18).  The 
terminal  descent  velocity  of  the  PRSLS  is  computed  from  the  well  known  parachute 
expression: 


^TERM~ 

N 


2  i^p-^Wpj^ 

TcrpC^p 


(27) 


where  Wp,  Wpis^r,  rp,  Cp,,  and  p  are  as  previously  defined  in  this  report. 

The  program  next  computes  the  initial  conditions  at  time  equals  zero,  i.e., 
when  the  valve  opens  and  retraction  begins.  Gas  expansion  work,  piston  velocity 
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and  the  displacements  of  the  piston,  payload,  and  parachute  confluence  point  are 
all  set  equal  to  zero  at  time  equals  zero.  Alfliough  the  payload  and  parachute 
confluence  point  are  both  initially  assumed  to  be  moving  downward  at  terminal 
descent  velocity,  as  an  acknowledgment  of  the  sign  convention  used  in  model 
development,  their  respective  velocities  are  taken  to  be  equal  to  plus  and  minus  the 
value  computed  using  equation  (27).  The  force  exerted  on  the  piston  by  the 
expanding  gas  is  computed  next  using  either  equation  (3)  or  equation  (5), 
depending  upon  wheflier  adiabatic  or  isothermal  gas  expansion  is  assumed  in  the 
analysis.  Using  the  model  inputs  supplied,  and  the  values  computed  for  the 
constant  terms,  the  force  acting  on  the  piston,  and  the  displacements  and  velocities 
of  the  three  masses,  the  program  is  now  able  to  compute  all  of  the  remaining 
forces,  displacements,  velocities  and  accelerations  at  time  zero  using  previously 
derived  expressions.  The  accelerations  of  the  parachute  confluence  point,  piston 
and  payload  are  first  computed  using  equations  (19),  (14a),  and  (14b)  respectively. 
With  the  motion  of  the  three  masses  now  completely  described  at  time  equals  zero, 
the  remaining  values  for  the  line  tension,  parachute  force,  retraction,  retraction  rate, 
and  retraction  acceleration  can  then  be  computed  respectively  using  equations  (12), 
(11),  (15a),  (15b),  and  (15c). 

Integration  of  the  PRSLS's  nonlinear  differential  equation  of  motion  in  time 
begins  by  first  computing  the  piston  force  and  gas  expansion  work  when  the  piston 
has  moved  a  distance.  Ay,  the  step  increment  chosen  by  the  program  user.  Force 
and  work  values  are  computed  using  either  equations  (3)  and  (4)  or  equations  (5) 
and  (6),  depending  upon  whether  adiabatic  or  isothermal  gas  expansion  has  been 
assumed  in  the  analysis.  Using  recurrence  formulas  analogous  to  those  shown  in 
equations  (20)  through  (26)  and  the  differential  equation  of  motion,  equation  (19), 
the  displacement,  velocity,  and  acceleration  of  the  parachute  confluence  point  are 
determined  iteratively  using  the  following  procedure. 

Assuming  a  value  for  the  change  in  time.  At,  equal  to  0.1  seconds,  A  and  B 
are  computed  at  the  beginning  time  station  of  flie  time  interval,  in  this  case,  time 
equals  zero,  using  equations  (22)  and  (25)  respectively.  One  will  note  that  all  of  the 
motion  values  needed  to  compute  A  and  B  are  known;  having  been  previously 
computed.  The  velocity,  displacement,  and  acceleration  of  the  parachute 
confluence  point,  at  the  end  of  the  time  interval,  are  then  computed  respectively 
using  equations  (26),  (24),  and  (19)  in  that  order.  Using  the  acceleration  value 
previously  computed,  updated  estimates  of  the  velocity,  displacement,  and 
acceleration  are  then  obtained  respectively  from  equations  (21),  (24),  and  (19).  The 
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latest  and  previously  computed  values  of  displacement  at  the  parachute  confluence 
point  are  then  compared.  If  the  absolute  value  of  the  difference  is  greater  than 
0.00001  inches,  velocity,  displacement  and  acceleration  values  are  updated,  once 
again  using  equations  (21),  (24),  and  (19),  and  confluence  point  displacements  are 
again  compared.  This  procedure  is  repeated  until  the  absolute  value  of  the 
difference  between  these  displacements  is  less  than  0.00001  inches.  At  that  point 
the  nonlinear  motion  equation  is  considered  solved  for  the  assumed  At. 


Once  the  differential  equation  of  motion  has  been  solved  for  the  assumed  At 
value,  the  piston  acceleration  and  piston  displacement  respectively  are  computed 
using  equation  (14a)  and  an  equation  analogous  to  equation  (23).  One  will  recall 
that  the  force  exerted  on  the  piston  by  the  expanding  gas,  F,  is  not  known  as  a 
function  of  time,  but  was  computed  from  previously  derived  thermodynamic 
relationships  at  a  specified  value  of  piston  displacement.  Therefore,  the  specified 
and  computed  values  of  the  piston  displacement  have  to  be  compared,  and  the 
change  in  time.  At,  adjusted  by  means  of  a  second  iterative  procedure  in  order  to 
reconcile  any  differences.  This  process  begins  with  an  examination  of  the 
difference  between  displacement  values,  i.e.,  computed  piston  displacement  minus 
the  specified  value.  If  this  value  is  greater  than  the  user-specified  step  increment. 
Ay,  the  At  value  assumed  in  the  solution  of  the  equation  of  motion  is  halved,  and 
the  equation  is  solved  once  again  using  the  iterative  procedure  previously 
described.  Specified  and  computed  values  of  the  piston  displacement  are  compared 
and  the  process  is  repeated  continuously  until  the  difference  is  less  than  the  step 
increment.  Ay.  At  that  point,  updated  values  of  At  are  obtained  by  means  of 
extrapolation  using  the  FORTRAN  formula: 


(  yoIFF 


I  Ay 


where  y^ipp  is  obtained  by  subtracting  the  specified  value  of  the  piston  displacement 
from  the  computed  value  during  each  iteration*.  The  equal  sign  in  the  above 
FORTRAN  formula  by  definition  means  "is  replaced  by".  Thus,  estimates  of  At 
are  either  increased  or  decreased  prior  to  again  solving  the  motion  equation, 
depending  respectively  upon  whether  one  has  undershot  (  ypipp  is  negative)  or 
overshot  (yoipp  is  positive)  the  value  of  At  during  the  previous  iteration. 


This  two-step  iterative  procedure  continues  until  the  absolute  value  of  y^iFF 
is  less  than  0.00001  inches.  At  that  point,  the  equation  of  motion  is  assumed  to  be 
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solved  for  that  time  step,  having  closely  coupled  the  kinetic,  kinematic  and 
thermodynamic  relationships  on  which  system  behavior  is  based.  The  remaining 
force,  displacement,  velocity,  and  acceleration  values  can  now  be  computed. 
Payload  acceleration  is  computed  first  using  equation  (14b).  Piston  and  payload 
velocity  are  then  computed  using  recurrence  formulas  analogous  to  equation  (20); 
followed  by  the  payload  displacement  using  a  recurrence  formula  similar  to 
equation  (23).  The  line  tension,  parachute  force,  retraction,  retraction  rate,  and 
retraction  acceleration  are  then  obtained  respectively,  once  again,  from  equations 
(12),  (11),  (15a),  (15b),  and  (15c).  Next,  tiie  elapsed  time  up  to  the  time  station 
is  computed  by  adding  At  to  the  elapsed  time  at  the  previous  time  station.  Finally, 
values  are  written  to  an  output  file,  the  piston  displacement  is  increased  by  the 
user-specified  step  increment.  Ay,  and  the  two-step  iterative  procedure,  described 
here,  is  carried  out  at  the  next  time  station,  using  the  displacements,  velocities,  and 
accelerations  from  the  most  recently  computed  time  station  as  initial  conditions. 
The  time  history  analysis  continues  one  step  at  a  time,  at  Ay  increments,  until 
either  the  time  reaches  the  maximum  input  value  specified  or  the  piston 
displacement  reaches  its  full  stroke.  The  analysis  is  terminated  once  one  of  these 
parameters  has  been  exceeded. 
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Comparison  of  Experimental  Data  with  Model  Predictions 

To  test  the  computer  model's  ability  to  predict  system  performance, 
FORTRAN  model  predictions  were  compared  with  experimental  data  obtained  from 
tests  conducted  of  a  prototype  PRSLS  at  Tustin  MCAS  in  April  1994.  A  complete 
verification  and  validation  of  the  computer  model  at  the  present  time,  however,  is 
not  possible.  This  is  due  to  the  fact  that  the  Tustin  experiments  remain  the  only 
instrumented  tests  conducted  to  date  on  a  PRSLS  of  this  kind,  and  the  amount  of 
data  gathered  was  limited.  Therefore,  only  a  cursory  check  of  the  computer 
model,  based  on  this  limited  test  data,  can  be  made  here. 

Data  recorded  in  the  20  April  1994  Tustin  experiment  are  presented  in  Figure 
3.  This  figure  shows  both  the  line  tension  and  payload  descent  velocity  plotted 
verses  time.  The  retraction  begins  at  approximately  7.5  seconds  into  the  system's 
descent  and  ends  at  approximately  9.0  seconds;  therefore,  the  total  retraction  takes 
approximately  1.5  seconds.  One  will  observe  that  at  7.5  seconds,  the  line  tension, 
or  tension  at  the  parachute  confluence  point,  rises  almost  vertically  from  62.5 
pounds,  the  suspended  weight,  to  about  100  pounds,  and  then  decreases  in  a  non¬ 
linear  manner  to  about  50  pounds  at  9.0  seconds.  The  payload  velocity, 
meanwhile,  decreases  nonlinearly  from  the  terminal  descent  velocity  of 
approximately  21  feet  per  second  at  the  beginning  of  retraction  to  around  one-third 
that  value  at  the  end. 

In  order  to  compare  FORTRAN  computer  model  predictions  with  the 
experimental  data  from  Tustin,  the  program  was  run  with  input  parameters  that 
closely  approximate  those  of  the  prototype  PRSLS  tested  at  Tustin^  Input 
parameters  specified  in  these  runs  were  as  follows: 

Initial  Pressure  -  220  psig. 

Tank  Volume  -  212  cubic  inches 
Piston  Radius  -  1.25  inches 
Specific  Heat  Ratio  -  1.4 

Weight  of  Piston,  Cross  Member,  and  Top  Pulleys  -  1.0  pound 
Weight  of  Payload  and  Retractor  -  62.5  pounds 
Piston  Stroke  -  24.0  inches 
Step  Increment  -  0.1  inches 

Gas  Expansion  Process:  (1)  Adiabatic  (2)  Isothermal 
Parachute  Diameter  -  14  feet 
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Parachute  Retractor  Performance  Test,  20Apr94 


Figure  3:  Tustin  Experimental  Data 
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Mechanical  Advantage  of  Pulley  System  -  10 
Drag  Coefficient  of  Parachute  -  0.75 
Maximum  Run  Time  -  4  seconds 

The  model  was  run  twice;  first  assuming  adiabatic  expansion  of  the  ideal  gas, 
and  then  isothermal  gas  expansion.  Although  adiabatic  expansion  of  the  gas  was 
expected  to  be  closer  to  actual  gas  behavior,  isothermal  gas  expansion  was  also 
investigated  to  determine  whether  appreciable  difference  between  the  two  processes 
can  be  expected  over  the  anticipated  short  duration  in  which  retraction  occurs. 
Model  results  were  summarized  in  a  series  of  plots  using  the  MATLAB  computer 
program^®.  A  fisting  of  this  MATLAB  program  is  shown  in  Appendix  B.  This 
program  allows  the  user  to  selectively  view  results  by  choosing  from  a  fist  of  20 
available  plots  of  various  model  output  data.  Model  predictions  of  payload  velocity 
and  fine  tension  vs.  time  are  presented  in  the  plots  shown  in  Figure  4.  Although 
no  attempts  were  made  to  refine  any  of  the  model  inputs,  there  appears  to  be  close 
correlation  between  computer  predictions  assuming  adiabatic  gas  expansion  and  the 
observed  test  behavior. 

Although  not  shown  in  the  fine  tension  plot,  one  will  observe  that  since  the 
system  is  initially  at  terminal  descent  velocity,  the  tension  at  the  parachute 
confluence  point  prior  to  retraction  is  equal  to  the  suspended  weight;  in  this  case, 
62.5  pounds.  Therefore,  at  time  equals  zero,  when  retraction  begins,  the  tension 
at  the  parachute  confluence  point  is  predicted  to  rise  vertically  from  the  suspended 
weight  up  to  approximately  108  pounds.  This  overprediction  of  the  peak  tension 
at  the  start  of  retraction  is  undoubtedly  due  to  the  instantaneous  start  assumed  in 
the  development  of  the  model.  The  predicted  tension  is  then  observed  to  decrease 
nonlinearly  to  54.6  pounds  at  the  end  of  the  retraction.  The  payload  velocity, 
meanwhile,  is  shown  to  decrease  nonlinearly  from  the  computed  terminal  descent 
velocity  of  256  inches  per  second  (21.3  feet  per  second)  to  a  minimum  of 
approximately  81.4  inches  per  second  (6.8  feet  per  second).  The  payload  velocity 
then  increases  slightly  before  the  piston  reaches  its  full  stroke.  Full  retraction  of 
the  parachute  is  predicted  to  take  place  in  just  under  1.5  seconds.  Thus,  close 
agreement  between  experimental  test  data  and  model  predictions  assuming  adiabatic 
gas  expansion  was  obtained  in  terms  of  not  only  the  magnitude  of  values,  but  also 
the  shapes  of  the  curves.  One  will  further  observe  that  significant  differences 
between  adiabatic  and  isothermal  computer  model  predictions  are  obtained  despite 
the  short  duration  of  the  retraction.  These  differences  become  more  pronounced 
toward  the  end  of  the  retraction;  a  product  of  the  isothermal  process's  larger  amount 
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of  available  energy  as  compared  with  that  of  the  adiabatic  process. 

The  reader’s  attention  is  directed  next  to  Figure  5,  showing  FORTRAN  model 
predictions  of  both  parachute  force  vs.  time,  and  payload  velocity  vs.  payload 
displacement.  Though  not  directly  used  in  the  comparisons  of  experimental  data 
with  model  predictions.  Figure  5  is  presented  here,  nonetheless,  by  way  of  example, 
to  provide  the  reader  a  greater  appreciation  of  the  model's  capabilities.  The  intent 
is  to  demonstrate  the  model's  usefulness  and  flexibility,  by  showing  how  the 
computer  model  can  be  used  by  system  developers  and  end  users  alike,  and  not  just 
as  a  research  tool.  Plots  of  the  parachute  force,  as  well  as  other  force  and 
acceleration  time  histories  provided  by  the  model,  for  example,  can  be  used  by 
system  developers  to  check  and  refine  the  design  of  the  retractor  mechanism  and 
parachutes,  prior  to  system  fabrication  and  testing.  Therefore,  much  of  the  costly, 
time  consuming  need  to  build,  instrument  and  test  prototypes,  could  be  eliminated 
during  the  Engineering  and  Manufacturing  Development  (EMD)  phase  of  system 
development,  by  using  the  model  as  a  design  tool.  Similarly,  given  a  PRSLS  of 
some  known  capacity,  one  can  quickly  compute  the  pressure  and  activation  height 
required  to  soft  land  a  payload  of  some  lesser  weight  from  plots  of  the  payload 
velocity  vs.  payload  displacement.  Therefore,  the  model  potentially  could  be  used 
as  a  simulation  tool  to  support  PRSLS  airdrop  operations  by  providing  the  end  user, 
i.e.,  troops  in  the  field,  Ae  ability  to  quickly  determine  soft  landing  requirements 
for  cargo  of  varied  weights. 
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Conclusions  and  Recommendations 


This  report  has  traced  the  development  of  a  predictive,  time  history  model 
of  the  NRDEC  Parachute  Retraction  Soft  Landing  System  (PRSLS);  describing  the 
model's  underlying  theory  and  current  status.  Based  on  comparisons  of 
experimental  data  with  model  predictions,  it  was  concluded  that  this  computer 
model,  in  its  current  form,  can  predict,  with  a  fair  degree  of  accuracy,  the  system's 
motion  and  forces  during  retraction.  The  model,  once  validated,  possesses  the 
flexibility  and  capability  to  serve  as  a  research,  design,  and  simulation  tool  that  will 
further  our  understanding  of  system  behavior,  simplify  the  design  process,  and  aid 
in  the  selection  of  pressures  and  activation  heights.  As  such,  this  model  represents 
an  enabling  technology  that  can  be  used  to  speed  introduction  of  the  PRSLS  into 
the  Army's  inventory.  Therefore,  it  is  recommended  that  more  experimental  data 
be  obtained  to  validate  the  model  over  a  range  of  payload  weights  and  pulley  ratios 
and/or  provide  a  foundation  for  further  model  refinements. 


This  document  reports  research  undertaken 
at  the  U.S,  Army  Soldier  Systems  Command, 
Natick  Research,  Development  and  Engineering 
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approved  for  publication. 
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Appendix  A 


c  PROGRAM :  ' RETRACTOR ' 

C 

Q  ***********************************************************^***^^^^ 

C  THIS  PROGRAM  MODELS  BOTH  THE  LANDING  RETRACTOR  MECHANISM  AND  THE 
C  PARACHUTE  CANOPY.  OUTPUT  INCLUDES  A  COMPLETE  DESCRIPTION  OF  THE 
C  DYNAMIC  MOTION  OF  THE  SYSTEM  INCLUDING  DISPLACEMENT,  VELOCITY  AND 
C  ACCELERATION  OF  THE  PISTON,  PAYLOAD  AND  CANOPY,  LINE  TENSION 
C  PARACHUTE  FORCE,  PISTON  WORK,  RETRACTION  RATE,  ETC.  VS.  TIME. 


★ 


IMPLICIT  DOUBLE  PRECISION (A-H, J-M, 0-Z) , INTEGER  (I,N) 

DIMENSION  Y(0:5000) ,YC (0:5000) ,YP (0:5000) ,VY (0:5000) ,VYC (0:5000 
DIMENSION  VYP(0:5000) ,AY(0:5000) ,AYC(0:5000)  ,AYP(0:5000) 
DIMENSION  F (0:5000) 

DIMENSION  PWORK(0 :5000) ,T(0:5000) , FPAR ( 0 : 5000) ,DL(0:5000) 
DIMENSION  VL (0:5000) ,AL (0:5000) , A (0:5000) ,B (0:5000) 


) 


TO  EXECUTE  PROGRAM  TYPE  'RETRACTOR'  AND  HIT  THE  RETURN  KEY 
TYPE  IN  VALUES  FOR  THE  FOLLOWING  VARIABLES  SEPARATED  EITHER  BY 
A  SPACE  OR  COMMA: 


INITIAL  PRESSURE  IN  PSI .  (GAUGE) 

INITIAL  VOLUME  IN  CUBIC  INCHES 
PISTON  RADIUS  IN  INCHES 
SPECIFIC  HEAT  RATIO,  k=Cp/Cv 

WEIGHT  OF  PISTON,  CROSS  MEMBER  AND  TOP  PULLEYS  IN  LBS. 
WEIGHT  OF  THE  PAYLOAD  AND  RETRACTOR  IN  LBS. 

PISTON  STROKE  IN  INCHES 
INCREMENT  OF  Y  IN  INCHES 

GAS  EXPANSION  PROCESS:  (1)  ADIABATIC  (2)  ISOTHERMAL 
PARACHUTE  DIAMETER  IN  FEET 


MECHANICAL  ADVANTAGE  OF  PULLEY  SYSTEM  (EVEN  INTEGER) 

DRAG  COEFFICIENT  OF  PARACHUTE 
MAXIMUM  RUN  TIME  IN  SECONDS 

******************************************************************* 


READ  *,PIN,VIN,R,k,WPIST,WP,  STROKE ,  DELTAY ,  IPROC,  DPAR,N,  CD,  TMAX 


WRITES  VARIABLE  VALUES  TO  OUTPUT  FILE 


50 


51 


52 


53 


54 


OPEN (24, FILE= ' OUTPUT ' ) 

WRITE (24, 50) PIN 

FORMAT ( ' INITIAL  PRESSURE  IN  PSI . ' / 

2  5X,F13.5) 

WRITE (24, 51) VIN 

FORMAT ( ' INITIAL  VOLUME  IN  CUBIC  INCHES ' / 

2  5X,F13.5) 

WRITE (24, 52 )R 

FORMAT ( ' PISTON  RADIUS  IN  INCHES ' / 

2  5X,F13.5) 

WRITE (24, 53 )k 

FORMAT (' SPECIFIC  HEAT  RATIO,  k=Cp/Cv'/ 

2  5X,F13.5) 

WRITE (24,54) WPIST 

FORMAT (' WEIGHT  OF  PISTON,  CROSS  MEMBER  AND  TOP  PULLEYS  IN  LBS  '/ 
2  5X,F13.5) 
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WRITE (24, 55 )WP 

55  FORMAT ( 'WEIGHT  OF  THE  PAYLOAD  AND  RETRACTOR  IN  LBS.  '/ 

2  5X,F13.5) 

WRITE (24,56) STROKE 

56  FORMAT (' PISTON  STROKE  IN  INCHES'/ 

2  5X,F13.5) 

WRITE (24,57) DELTAY 

57  FORMAT  (  '  INCREMENT  OF  Y  IN  INCHES '  / 

2  5X,F13.5) 

IF(IPROC.EQ.l)  THEN 
WRITE (24, 58) 

58  FORMAT ( 'ADIABATIC  GAS  EXPANSION  PROCESS') 

GO  TO  61 

ELSE 

CONTINUE 

ENDIF 

IF(IPROC.EQ.2)  THEN 
WRITE (24, 60) 

60  FORMAT  (' ISOTHERMAL  GAS  EXPANSION  PROCESS') 

GO  TO  61 

ELSE 

GO  TO  31 
ENDIF 

61  WRITE(24, 62)DPAR 

62  FORMAT ( 'PARACHUTE  DIAMETER  IN  FEET'/ 

2  5X,F13.5) 

WRITE (24, 63 )N 

63  FORMAT ( 'MECHANICAL  ADVANTAGE  OF  PULLEY  SYSTEM'/ 

2  5X,I4) 

WRITE (24, 64) CD 

64  FORMAT ('DRAG  COEFFICIENT  OF  PARACHUTE'/ 

2  5X,F13.5) 

WRITE(24,65)TMAX 

65  FORMAT (' MAXIMUM  RUN  TIME  IN  SECONDS'/ 

2  5X,F13.5) 

DEFINE  AND  COMPUTE  CONSTANT  TERMS  USED  IN  THE  DIFFERENTIAL 
EQUATIONS  OF  MOTION 

g=32. 174*12.0 
pi=3 .141592654 
PA=14 . 7 
RPAR=DPAR/2 . 

WAIR=0. 07648 
rho=WAIR/32 .174 

VTERM=12 . 0*SQRT (2 . 0*WP/ (pi*RPAR**2 . *CD*rho) ) 

WP=WP-WPIST 

MP=WP/g 

MPIST=WPIST/g 

MC= (2 . *pi*rho*RPAR**3 . ) / (3 . *12 . ) 

Cl=pi*RPAR**2.*CD*rho/ (2.*144. ) 

C2=N**2*C1*MP+ (N+1) *C1*MPIST 

C3= (N+1) *MC*MPIST+MP*MPIST+N**2*MC*MP 

C4= (N-1) *g*MPIST*MP 

C5=N*MP+MPIST 

PIN=PIN+PA 
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BEGINNING  OF  COMPUTATIONS  AT  FIRST  TIME  STATION,  Y{0)=0.,  TIME=0 . 

YPIST=0 . 

1=0 

TIME=0 . 

Y(I)=YPIST 
YC(I)=0. 

YP(I)=0. 

VY{I)=0. 

VYC{I) =VTERM 
VYP(I)=-VTERM 
PWORK(I)=0. 

IF(IPROC.EQ.l)  THEN 

F (I) = (PIN*VIN**k*pi*R**2 . / (VIN+pi*R**2 . *Y(I) ) **k) -PA*pi*R**2 . 
GO  TO  11 
ELSE 
GO  TO  10 
ENUIF 

0  F(I)=(PIN*VIN*pi*R**2./ (VIN+pi*R**2.*Y(I) ) ) -PA*pi*R**2 . 

1  AYC(I) =(C5*F(I)-C2*VYC(I) **2.-C4) /C3 

AY ( I )  =  { F ( I ) -N*C1*VYC ( I ) *  *2 . -N*MC*AYC ( I )  -WPIST)  /MPIST 
AYP(I) =( (N+l) *C1*VYC(I) **2.+(N+l) *MC*AYC (I) -F (I) -WP) /MP 
T(I)=C1*VYC{I)**2.+MC*AYC(I) 

FPAR(I)=C1*VYC(I) **2 . 

DL(I)=YC{I)+YP(I) 

VL(I)=VYC(I)+VYP{I) 

AL(I)=AYC(I)+AYP(I) 

GO  TO  30 

END  OF  COMPUTATIONS  AT  FIRST  TIME  STATION,  Y(0)=0.,  TIME=0 . 

20  IF(IPROC.EQ.l)  THEN 

F(I)=(PIN*VIN**k*pi*R**2 . / (VIN+pi*R**2 . *Y(I) ) **k) -PA*pi*R**2 . 
PWORK(I) = (PIN*VIN**k* (VIN+pi*R**2 . *Y(I) ) ** (l-k) -PIN*VIN) / (1-k) 
PWORK ( I ) =PWORK ( I ) -PA*pi*R*  *2 . *Y ( I ) 

GO  TO  22 
ELSE 

GO  TO  21 
ENDIF 

21  F (I) = (PIN*VIN*pi*R**2 . / (VIN+pi*R**2 . *Y(I) ) ) -PA*pi*R**2 . 

PWORK  (I)  =PIN*VIN*DLOG  (  (VIN+pi*R**2  .  *Y{I)  )  /VIN) 

PWORK ( I ) =PWORK ( I ) -PA*pi *R*  *2 . *Y ( I ) 

22  DELTAT=0 . 1 

25  A(I-1)=VYC(I-1)+AYC(I-1) *DELTAT/2. 

B ( I-l ) =YC ( I-l ) +2 . *VYC ( I-l ) *DELTAT/3 . +AYC ( I-l ) *DELTAT*  *2 . / 6 . 

VYC  ( I )  =VYC  ( I-l )  +AYC  ( I-l )  *DELTAT 
YC{I)=B(I-1)+VYC(I) *DELTAT/3 . 
AYC(I)=(C5*F(I)-C2*VYC(I)**2.-C4)/C3 
PREVYC=YC(I) 

23  VYC(I)=A(I-l)+AYC(I)*DELTAT/2. 

YC ( I ) =B ( I-l ) +VYC ( I ) *DELTAT/3 . 

AYC{I)=(C5*F(I)-C2*VYC{I) **2.-C4)/C3 

IF( (DABS(YC(I)-PREVYC) ) .LT. 0.00001)  THEN 

GO  TO  24 

ELSE 
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PREVYC=YC(I) 

GO  TO  23 
ENDIF 

24  AY{I)=(F(I) -N*C1*VYC(I) **2 . -N*MC*AYC (I) -WPIST) /MPIST 

Y ( I ) =Y ( I-l ) +VY ( I-l ) *DELTAT+ ( 2 . *AY { I-l ) +AY ( I ) ) *DELTAT*  *2 . / 6 . 

YDIFF=Y(I)-YPIST 

IF(YDIFF  .GE.  DELTAY)  THEN 

DELTAT=DELTAT/2 . 

GO  TO  25 
ELSE 

CONTINUE 

ENDIF 

IF( (YDIFF.lt. DELTAY) .AND. ( {DABS(YDIFF) ) .GT. 0.00001) )  THEN 
DELTAT=DELTAT-DELTAT* ( YDIFF /DELTAY ) 

GO  TO  25 
ELSE 

VY ( I ) =VY ( I-l ) +DELTAT/2 . * (AY ( I ) +AY ( I-l ) ) 

AYP(I)=( (N+1)*C1*VYC(I) **2.+(N+l) *MC*AYC(I) -F(I) -WP) /MP 
VYP(I)=VYP(I-1)+(AYP(I-1)+AYP(I) ) *DELTAT/2 . 

YP(I)=YP(I-1)+VYP(I-1) *DELTAT+(2.*AYP(I-1)+AYP(I) ) *DELTAT**2 . /6 . 
T(I) =C1*VYC(I) **2 .+MC*AYC(I) 

FPAR(I)=C1*VYC(I) **2 . 

DL(I)=YC(I)+YP(I) 

VL  (I)  =VYC  (I) +VYP  (I) 

AL(I)=AYC(I)+AYP(I) 

TIME=TIME+DELTAT 

ENDIF 

WRITE  RESULTS  TO  OUTPUT  FILE  AND  COMPARE  TIME  AND  PISTON 
DISPLACEMENT  RESPECTIVELY  TO  MAXIMUM  TIME  AND  PISTON  STROKE 

30  IF (TIME  .GT.  TMAX)  GO  TO  31 

WRITE(24,66)I,F(I) , DELTAT, TIME, Y ( I) ,YC(I) ,YP(I) ,AY(I) ,AYC(I) , 
2AYP(I) ,VY(I) ,VYC(I) ,VYP(I) ,PWORK(I) ,T(I) ,FPAR(I) , 

3DL(I) ,VL(I) ,AL(I) 

66  FORMAT(2X,I5,18 (2X,E13 .5) ) 

IF(Y(I)  .LT.  STROKE)  THEN 
1=1+1 

YPIST=YPIST+DELTAY 
Y(I)=YPIST 
GO  TO  20 
ELSE 

CONTINUE 

ENDIF 

31  CLOSE (24 , STATUS= ' KEEP ' ) 

STOP 

END 
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MATUVB  Program  Iiis'ting 


clear 

load  -ascii  PLRl.mat 
load  -ascii  PLR2.mat 
who 
clg 

hold  off 

k  =  menu ( ' CHOOSE  PLOT ’ , . . . 

' PISTON  ACCELERATION  VS .  TIME ' , . . . 

' CONFLUENCE  POINT  ACCELERATION  VS .  TIME ' , . . . 

' PAYLOAD  ACCELERATION  VS .  TIME ' , . . . 

' PISTON  VELOCITY  VS .  TIME ' , . . . 

' CONFLUENCE  POINT  VELOCITY  VS .  TIME ' , . . . 

' PAYLOAD  VELOCITY  VS .  TIME ' , . . . 

' PISTON  DISPLACEMENT  VS .  TIME ' , . . . 

' CONFLUENCE  POINT  DISPLACEMENT  VS .  TIME ' , . . . 

’ PAYLOAD  DISPLACEMENT  VS ,  TIME ’ , . . . 

■PISTON  WORK  VS.  RETRACTION' , . . . 

' PISTON  FORCE  VS .  DISPLACEMENT ' , . . . 

■PISTON  WORK  VS.  DISPLACEMENT' , . . . 

■  PISTON  FORCE  VS .  TIME ' ,  . . . 

■  PISTON  WORK  VS .  TIME ' ,  . . . 

■  CONFLUENCE  POINT  TENSION  VS .  TIME ' ,  . . . 

■PARACHUTE  FORCE  VS.  TIME' , . . . 

■  RETRACTION  VS .  TIME ' ,  . . . 

■  RETRACTION  VELOCITY  VS .  TIME ' ,  . . . 

■  RETRACTION  ACCELERATION  VS .  TIME  '  ,  .  .  . 

■  PAYLOAD  VELOCITY  VS .  PAYLOAD  DISPLACEMENT ' ,  .  . . 

■  QUIT ■ ) 
if  k  ==  1 
clg 

plot (PLRl ( : , 4) , PLRl ( : , 8) , ' -r ' , PLR2 ( : , 4) , PLR2 ( : , 8 ) , ' -w' ) 
title ( 'PISTON  ACCELERATION  VS.  TIME') 
xlabel ( ■ TIME  SECONDS ' ) 

ylabel ( 'ACCELERATION  INCHES  PER  SEC  2') 

% 

elseif  k  ==  2 
clg 

plot (PLRl ( : , 4) , PLRl ( : , 9) ,  ' -r ' , PLR2 { : , 4)  ,  PLR2 ( :  ,  9)  ,  ' -w' ) 
title ( 'CONFLUENCE  POINT  ACCELERATION  VS.  TIME') 
xlabel ( ' TIME  SECONDS ' ) 

ylabel ( 'ACCELERATION  INCHES  PER  SEC  2') 

% 

elseif  k  ==  3 
clg 

plot (PLRl ( : ,4) ,PLR1 ( : , 10) , ' -r ' , PLR2 ( : ,4) , PLR2 ( : , 10) , ' -w' ) 
title ( ’ PAYLOAD  ACCELERATION  VS .  TIME ' ) 
xlabel ( ' TIME  SECONDS ' ) 

ylabel ( 'ACCELERATION  INCHES  PER  SEC  2') 

% 

elseif  k  ==  4 
clg 

plot (PLRl ( : , 4) , PLRl ( : , 11) , ' -r ' , PLR2 ( : , 4) , PLR2 ( : , 11) , ’ -w' ) 
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title ( 'PISTON  VELOCITY  VS.  TIME') 

xlabel ( ' TIME  SECONDS ' ) 

ylabel  ( 'VELOCITY  INCHES  PER  SEC) 

% 

elseif  k  ==  5 
clg 

plot (PLRl ( : , 4 ) , PLRl ( : , 12 ) , ' -r ' , PLR2 ( : , 4) , PLR2 ( : , 12 ) , ' -w' ) 

title ( 'CONFLUENCE  POINT  VELOCITY  VS.  TIME') 

xlabel ( ' TIME  SECONDS ' ) 

ylabel  ( 'VELOCITY  INCHES  PER  SEC) 

% 

elseif  k  ==  6 
clg 

plot (PLRl ( : , 4 ) , PLRl { : , 13 ) , ' -r ' , PLR2 ( : , 4 ) , PLR2 ( : , 13 ) , ' -W ) 

title ( 'PAYLOAD  VELOCITY  VS.  TIME') 

xlabel ( ' TIME  SECONDS ' ) 

ylabel ( 'VELOCITY  INCHES  PER  SEC ' ) 

% 

elseif  k  ==  7 
clg 

plot (PLRl ( : ,4) ,PLR1 ( : ,5) , ' -r' , PLR2 { : ,4) , PLR2 { : ,5) , ' -W ) 
title ( 'PISTON  DISPLACEMENT  VS.  TIME') 
xlabel ( ' TIME  SECONDS ' ) 
ylabel ( ' DISPLACEMENT  INCHES ' ) 

% 

elseif  k  ==  8 
clg 

plot (PLRl ( : ,4) ,PLR1( : , 6) , ' -r ' , PLR2 ( : ,4) ,PLR2 ( : , 6) , '-W ) 
title ( 'CONFLUENCE  POINT  DISPLACEMENT  VS.  TIME') 
xlabel ( ' TIME  SECONDS ' ) 
ylabel { ' DISPLACEMENT  INCHES ' ) 

% 

elseif  k  ==  9 
clg 

plot (PLRl ( : ,4) ,PLR1( : ,7) , '-r' , PLR2 ( : ,4) , PLR2 ( : ,7) , '-W ) 
title ( 'PAYLOAD  DISPLACEMENT  VS.  TIME') 
xlabel ( ' TIME  SECONDS ' ) 
ylabel ( ' DISPLACEMENT  INCHES ' ) 

% 

elseif  k  ==  10 
clg 

plot (PLRl ( : ,17) ,PLR1( : ,14) , ' -r' ,PLR2 ( : ,17) , PLR2 ( : ,14) , '-w' 
title ( 'PISTON  WORK  VS.  RETRACTION') 
xlabel ( ' RETRACTION  INCHES ' ) 
ylabel ( ' WORK  INCH | POUNDS ' ) 

% 

elseif  k  ==  11 
clg 

plot(PLRl( : ,5) ,PLR1( : ,2) , '-r' , PLR2 ( : ,5) , PLR2 ( : ,2) , '-w' ) 
title ( 'PISTON  FORCE  VS.  DISPLACEMENT') 
xlabel ( ' DISPLACEMENT  INCHES ' ) 
ylabel ( ' FORCE  POUNDS ' ) 

% 

elseif  k  ==  12 
clg 

plot (PLRl ( : , 5 ) , PLRl ( : , 14 ) ,  ' -r ' , PLR2 ( : , 5 ) , PLR2 ( : , 14 ) ,  '  -w'  ) 
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title (' PISTON  WORK  vs .  DISPLACEMENT') 
xlabel { ■ DISPLACEMENT  INCHES ' ) 
ylabel { ' WORK  INCH | POUNDS ' ) 

% 

elseif  k  ==  13 
clg 

plot { PLRl ( : , 4 ) , PLRl ( : , 2 ) , ' -r ’ , PLR2 ( : , 4 ) , PLR2 ( : , 2 )  ,  ' -w ' ) 
title { ' PISTON  FORCE  VS .  TIME ' ) 
xlabel ( ■ TIME  SECONDS ' ) 
ylabel ( ' FORCE  POUNDS ' ) 


elseif  k  ==  14 
clg 

plot (PLRl ( : , 4) , PLRl ( : , 14) , ’ -r ' , PLR2 ( : , 4) , PLR2 ( : , 14) , ' -w' ) 
title (' PISTON  WORK  VS.  TIME') 
xlabel ( ' TIME  SECONDS ' ) 
ylabel ( ' WORK  INCH | POUNDS ' ) 

% 

elseif  k  ==  15 
clg 

plot (PLRl ( : , 4) , PLRl ( : , 15) , ' -r ' , PLR2 ( : , 4) , PLR2 ( : , 15) , ' -w' ) 
title ( ' CONFLUENCE  POINT  TENSION  VS .  TIME ' ) 
xlabel ('TIME  SECONDS') 
ylabel ( ' TENSION  POUNDS ' ) 

% 

elseif  k  ==  16 
clg 

plot (PLRl (: ,4) ,PLR1 ( : ,16) ,  '-r' , PLR2 ( : ,4) , PLR2 ( : , 16) ,  '-W  ) 
title (' PARACHUTE  FORCE  VS.  TIME') 
xlabel ( ' TIME  SECONDS ' ) 
ylabel ( ' FORCE  POUNDS ' ) 

% 

elseif  k  ==  17 
clg 

plot (PLRl ( : , 4) , PLRl ( : , 17) , ’ -r ' , PLR2 ( : , 4) , PLR2 ( : , 17) , ' -w' ) 
title ( ' RETRACTION  VS .  TIME ' ) 
xlabel ( ' TIME  SECONDS ' ) 
ylabel ( ' RETRACTION  INCHES ' ) 

% 

elseif  k  ==  18 
clg 

plot (PLR1( : ,4) ,PLR1 ( : ,18) , '-r' , PLR2 ( : ,4) , PLR2 ( : , 18) , '-w' ) 

title ( 'RETRACTION  VELOCITY  VS.  TIME') 

xlabel ( ' TIME  SECOND ' ) 

ylabel ( ' VELOCITY  INCHES  PER  SEC ' ) 

% 

elseif  k  ==  19 
clg 

plot (PLR1( : ,4) ,PLR1 ( : ,19) ,  '-r' , PLR2 ( : ,4) , PLR2 ( : ,19) ,  '-w'  ) 
title ( 'RETRACTION  ACCELERATION  VS.  TIME') 
xlabel ('TIME  SECONDS') 

ylabel ( ' ACCELERATION  INCHES  PER  SEC  2 ' ) 

% 

elseif  k  ==  20 
clg 

plot (PLRl ( ; , 7) , PLRl ( : , 13) ,  ' -r ' , PLR2 ( : , 7) , PLR2 ( : , 13 ) ,  ' -w'  ) 
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title {' PAYLOAD  VELOCITY  VS.  PAYLOAD  DISPLACEMENT') 
xlabel ( ' DISPLACEMENT  INCHES ' ) 
ylabeK 'VELOCITY  INCHES  PER  SEC) 
end 
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List  of  Symbols 


Cp  parachute  drag  coefficient 

F  force  exerted  on  piston  by  expanding  gas 

Fparachute  parachute  force 

H  enthalpy 

L  retraction,  i.e.,  L=Ny=yc+yp 

Me  parachute  apparent  mass 

Mp  mass  of  payload  and  retractor  components  that  move  with  payload 

MpisT  mass  of  piston  and  retractor  components  that  move  with  piston 
N  mechanical  advantage 

P  pressure 

Pj  initial  tank  pressure 

Pj  initial  pressure 

P2  final  pressure 

Q  heat 

R  universal  gas  constant,  where  R  =  Cp  -  c^ 

S  piston  stroke 

T  temperature 

Tj  initial  temperature 

T2  final  temperature 

Tc  tension  at  parachute  confluence  point 

U  internal  energy 

V  volume 

Vj  initial  volume 

V2  final  volume 

Vj  tank  volume 

W  work 

Wq  adiabatic  work 

isothermal  work 

Wq  weight  of  air  in  the  parachute;  assumed  negligible 

Wp  weight  of  payload  and  retractor  components  that  move  with  payload 

WpisT  weight  of  piston  and  retractor  components  that  move  with  piston 
Cp  specific  heat  at  constant  pressure 

c„  specific  heat  at  constant  volume 

k  specific  heat  ratio,  Cp/c^ 

n  number  of  moles  of  gas 

r  piston  radius 


t 

'^TERM 

X 

y 

Yc 

yp 

Ydiff 
AH 


AU 

At 

Ay 

P 


ideal  Gas 
Ideal  Gas 


nominal  parachute  radius 
time 

terminal  descent  velocity 
some  generalized  displacement 
piston  displacement 

displacement  at  parachute  confluence  point 
payload  displacement 

computed  minus  specified  piston  displacement 

change  in  enthalpy  of  an  ideal  gas 

change  in  internal  energy  of  an  ideal  gas 

change  in  time  over  interval 

input  step  increment  of  piston  displacement 

density  of  air 


Throughout  this  report,  a  dot  or  double  dot  notation  over  any  of  the 
displacements  (x,  y,  y^  ,  yp  or  L)  will  be  used  to  designate  the  first  and  second 
derivatives  of  that  displacement  with  respect  to  time,  or  in  other  words,  the  velocity 
and  acceleration. 
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